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ABSTRACT 

The cumulative comoving number-density of galaxies as a function of stellar mass or 
central velocity dispersion is commonly used to link galaxy populations across dif¬ 
ferent epochs. By assuming that galaxies preserve their number-density in time, one 
can infer the evolution of their properties, such as masses, sizes, and morphologies. 
However, this assumption does not hold in the presence of galaxy mergers or when 
rank ordering is broken owing to variable stellar growth rates. We present an analysis 
of the evolving comoving number density of galaxy populations found in the Illustris 
cosmological hydrodynamical simulation focused on the redshift range 0 < z < 3. 
Our primary results are as follows: 1) The inferred average stellar mass evolution 
obtained via a constant comoving number density assumption is systematically bi¬ 
ased compared to the merger tree results at the factor of ~2(4) level when tracking 
galaxies from redshift z = 0 out to redshift z = 2(3); 2) The median number density 
evolution for galaxy populations tracked forward in time is shallower than for galaxy 
populations tracked backward in time; 3) A similar evolution in the median number 
density of tracked galaxy populations is found regardless of whether number density 
is assigned via stellar mass, stellar velocity dispersion, or dark matter halo mass; 4) 
Explicit tracking reveals a large diversity in galaxies’ assembly histories that cannot 
be captured by constant number-density analyses; 5) The significant scatter in galaxy 
linking methods is only marginally reduced by considering a number of additional 
physical and observable galaxy properties as realized in our simulation. We provide 
fits for the forward and backward median evolution in stellar mass and number den¬ 
sity for use with observational data and discuss the implications of our analysis for 
interpreting multi-epoch galaxy property observations as related to galaxy evolution. 

Key words: methods: numerical - cosmology: theory - cosmology: galaxy formation 
- galaxies: abundances 


Modeling galaxy evolution based on observational data re¬ 
quires a method for linking galaxy populations between 
different epochs. Establishing direct progenitor-descendant 
links would allow for reconstruction of the mass, size, star 
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1 INTRODUCTION 


formation rate, color, and morphology evolution (among 
other things) of galaxies directly from observational data. 
However, formulating a method that accurately links pro¬ 
genitor and descendant galaxy populations is non-trivial. 
Incorrectly linking observed galaxy populations between dif¬ 
ferent epochs results in errors in the inferred evolutionary 
tracks. This effect has come to be known as “progenitor 
bias” (e.g., van Dokkum & Franx 1996; Saglia et al. 2010) 
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and is a well-known effect that needs to be addressed in 
order to infer the mass (van Dokkum et al. 2010; Bram- 
mer et al. 2011; Patel et al. 2013; van Dokkum et al. 2013), 
size (Fan et al. 2008; Valentinuzzi et al. 2010; Carollo et al. 
2013; Patel et al. 2013; van Dokkum et al. 2013; Morishita 
et al. 2015), star formation rate (Papovich et al. 2011) and 
morphology evolution (van Dokkum & Franx 2001; Daddi 
et al. 2005). 

Several approximations have been employed to estimate 
progenitor-descendant linking of galaxy populations in order 
to minimize or remove progenitor bias. One method that has 
been applied is to select rare or distinct galaxy populations 
which one might reasonably be able to recover at differ¬ 
ent observational epochs. Such analysis has been applied to 
examine the redshift-dependent properties of brightest clus¬ 
ter galaxies (BCGs; e.g., Butcher & Oemler 1984; Aragon- 
Salamanca et al. 1993; Lidman et al. 2012; Vulcani et al. 
2014) as well as to populations of massive early-type galax¬ 
ies (ETGs; e.g., Daddi et al. 2005; Trujillo et al. 2007; Barro 
et al. 2013). The premise behind this linking metric is that 
both massive ETGs and BCGs observed at high redshift will 
remain massive ETGs and BCGs into the low-redshift uni¬ 
verse. This assumption is generally true. However, caution 
must be taken as the fraction of galaxies that are massive 
and quenched grows with time. This will drive an increase 
in the number of massive ETGs and BCGs that exist in the 
local universe compared to the high-redshift universe, and 
therefore create contamination in the inferred progenitor- 
descendant galaxy populations. It has been argued in Car¬ 
ollo et al. (2013) that this effect alone can result in the in¬ 
ferred size evolution of massive ETGs (but see also Belli 
et al. 2014; Keating et al. 2015, for careful analysis that 
indicates progenitor bias is insufficient to fully explain the 
observed size evolution). The current uncertainty that sur¬ 
rounds the size evolution of ETGs is dominated by a lack 
of theoretical understanding of how to link high and low 
redshift galaxy populations. 

Another approach - which is the focus of this paper - 
is to assume that progenitor and descendant galaxy popu¬ 
lations can be linked based on their cumulative comoving 
number density (e.g., Wake et al. 2006; van Dokkum et al. 
2010; Papovich et al. 2011; Brammer et al. 2011; Patel et al. 
2013; van Dokkum et al. 2013). This method is appealing be¬ 
cause it provides a straightforward way to infer information 
about galaxy evolution directly from multi-epoch cumula¬ 
tive stellar mass functions. Because no assumptions about 
galactic characteristics (e.g., massive and quenched) are re¬ 
quired, this method has been employed to infer the mass 
evolution of Milky Way progenitors (Patel et al. 2013; van 
Dokkum et al. 2013), the evolution of red/quenched galaxy 
fractions (Brammer et al. 2011), and the size and morpholog¬ 
ical evolution of massive galaxies (van Dokkum et al. 2010). 

Linking galaxy populations in this fashion implicitly in¬ 
volves two assumptions: (i) that the total number (density) 
of galaxies is conserved and (ii) that galaxies maintain their 
rank order. If both of these assumptions are true, then the 
most(least) massive galaxies at some initial redshift (e.g., 
2 = 2) will still be the most (least) massive galaxies at any 
other redshift (e.g., 2 = 0). The two primary issues with this 
approach are that galaxies undergo merger events which can 
reduce their number density, and that individual galaxies 
may have scattered/stochastic growth histories which can 


break their rank ordering. Mergers will drive an evolution in 
the number density of any given galaxy population by simply 
changing the number density of the total galaxy population 
with time, even if galaxies remain well rank-ordered. On the 
other hand, scattered growth rates will cause a population 
of galaxies with similar initial stellar mass at some initial 
redshift (e.g., 2 = 2) to have some potentially wide distri¬ 
bution of masses by the time they evolve to some new ob¬ 
servational epoch (e.g., 2 = 0). To avoid unnecessary biases 
when applying constant comoving number-density analysis 
to observational datasets, it is important to understand the 
extent to which galaxy mergers or variable galaxy growth 
rates impact number-density matching. 

Capturing galaxy mass assembly, including both merg¬ 
ers and variable mass growth rates, can be done in detail by 
employing numerical galaxy formation simulations. The ef¬ 
ficiency of constant comoving number-density selections has 
been studied using semi-analytic models (Leja et al. 2013; 
Mundy et al. 2015) and with abundance matching (Behroozi 
et al. 2013). Leja et al. (2013) employed the Guo et al. 
(2011) semi-analytic models based on the Millennium Simu¬ 
lation (Springel et al. 2005b) to compare tracked mass evo¬ 
lution against assumed mass evolution based on number- 
density selections. They found that a constant comoving 
number-density selection yielded inferred median descen¬ 
dant masses of high-redshift galaxy populations which dif¬ 
fered by 40% from the actual descendant masses. By apply¬ 
ing a correction to account for the scatter in galaxy growth 
rates and mergers they were able to reduce the mass offset 
error from a number-density selection to 12%. However, even 
with such a correction, significant scatter remains among the 
inferred growth rates. 

Mundy et al. (2015) used several different semi-analytic 
models also based on the Millennium Simulation (Springel 
et al. 2005b) to provide fitting functions that describe 
the “recovery fraction” of galaxies with time. They found 
that for their lowest number-density bin (corresponding to 
the most massive/rare haloes) a constant number-density 
selection recovered roughly 30% of available descendants, 
with the majority of the recovered galaxy population be¬ 
ing contamination. Behroozi et al. (2013) used an abun¬ 
dance matching based on the Bolshoi Simulation (Klypin 
et al. 2011) to infer the number-density evolution of galax¬ 
ies. They found a median evolution in the number density 
for all galaxy mass bins, and provided fitting functions for 
the median evolution. These fitting functions have been ap¬ 
plied in Marchesini et al. (2014), where the mass and star 
formation rate histories of massive galaxies were inferred us¬ 
ing the prescribed number density evolution tracks. 

Our work builds on these previous theoretical studies 
by examining the constant comoving number-density selec¬ 
tion method using a full volume hydrodynamical simula¬ 
tion. Multi-epoch constant comoving number-density selec¬ 
tion method has only recently been examined using hydro- 
dynamical simulations at high redshifts (Jaacks et al. 2015), 
with limited study at 2 < 3. To the extent that both hy¬ 
drodynamical simulations and semi-analytic methods accu¬ 
rately reproduce the evolution of the galaxy stellar mass 
function (see Somerville & Dave 2014, for a comparison and 
discussion), it is unlikely that full hydrodynamical simula¬ 
tions would yield very different results compared to their 
semi-analytic model counterparts. However, there are some 
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subtle but important issues that we might expect to yield 
concrete differences between the previous analyses and what 
we present in this paper. 

Specifically, the cosmology used in the original Millen¬ 
nium Simulation (Springel et al. 2005b) was based on the 
WMAP-1 data release (Bennett et al. 2003), and therefore 
applied a ag value that is higher than what is currently 
accepted. This has the consequence of boosting the halo- 
halo merger rate, which is of direct relevance to galaxy 
number-density analysis (Leja et al. 2013; Mundy et al. 
2015). A direct comparison of the halo-halo merger rates 
between the Millennium Simulation and Illustris shows that 
the merger rates are in good agreement despite differences 
in the adopted cosmological parameters (Rodriguez-Gomez 
et al. 2015). The galaxy-galaxy merger rates used in semi- 
analytic models can be significantly offset from the galaxy- 
galaxy merger rates found in Illustris (Rodriguez-Gomez 
et al. 2015) which could impose differences on the number 
density evolution. However, given that the results presented 
in Behroozi et al. (2013) (which used updated cosmological 
parameters) generally agreed with Leja et al. (2013), it is 
unlikely that updated cosmological parameters alone would 
drive major differences in the number-density evolution of 
galaxies. 

Perhaps more importantly, while semi-analytic models 
and hydrodynamical simulations both attempt to include a 
similar array of physical processes (gas cooling, star forma¬ 
tion, feedback, etc.), in detail the growth rates of galaxies are 
calculated in very different ways. It is therefore not guaran¬ 
teed that the scattered growth rates that can contribute to 
breaking galaxy rank order will be equally realized in semi- 
analytic models and hydrodynamical simulations. Moreover, 
while semi-analytic models are able to infer a wide range of 
internal galactic properties (e.g., galaxy size, bulge-to-disk 
ratio, etc.), all of these properties are managed and evolved 
at a sub-grid level. Although the hydrodynamical simula¬ 
tions employed in this paper also apply sub-grid models to 
manage many aspects of galaxy formation physics (e.g., star 
formation, ISM gas phase structure, etc.) some important 
galactic characteristics such as stellar velocity dispersion 
can be self-consistently evolved. Finally, the general anal¬ 
ysis presented in this paper complements nicely the detailed 
studies of specific galaxy populations that form in our simu¬ 
lations and have been compared against observations, such 
as the formation and evolution of compact massive galax¬ 
ies (Wellons et al. 2015a,b). 

The paper is structured as follows: In Section 2 we dis¬ 
cuss our methods including a description of the simulations 
that have been employed, the construction of galaxy prop¬ 
erty catalogs, and the merger trees that form the core of our 
analysis. In Section 3 we present the simulated multi-epoch 
cumulative stellar mass function and the corresponding in¬ 
ferred stellar mass evolution. We compare this against the 
actual tracked stellar mass evolution found in the simula¬ 
tion. We introduce useful fitting formulae here, which de¬ 
scribe the mass and number-density evolution for the explic¬ 
itly tracked galaxy population. In Section 4 we present the 
simulated multi-epoch cumulative velocity dispersion func¬ 
tion and consider whether velocity dispersion might act as 
a better galactic parameter for linking galaxy populations. 
In Section 5 we discuss our results, with a focus on under¬ 
standing the empirical origin of the trends we find in our 


simulations and exploring implications for using number- 
density analyses on observational data sets. We conclude in 
Section 6. 


2 METHODS 

In this paper we use the Illustris simulation to study the rela¬ 
tionship between the stellar mass, dark matter mass, central 
stellar velocity dispersion, and number-density evolution of 
galaxies. Full details of the Illustris project can be found 
in Vogelsberger et al. (2014a,b); Genel et al. (2014). 

Briefly, the Illustris simulation is a cosmological hydro¬ 
dynamical simulation run in a periodic box of size L = 
106.5 Mpc. Illustris was run using the AREPO simulation 
code (Springel 2010) using a physical setup that includes 
gravity, hydrodynamics, radiative cooling of gas (Katz et al. 
1996), star formation with associated feedback (Springel & 
Hernquist 2003), mass and metal return to the interstel¬ 
lar medium from aging stellar populations (Wiersma et al. 
2009; Vogelsberger et al. 2013), and supermassive black hole 
growth with associated feedback (Di Matteo et al. 2005; 
Springel et al. 2005a; Sijacki et al. 2007; Vogelsberger et al. 
2013; Sijacki et al. 2014). The feedback models employed 
in the simulation were chosen to match the redshift z = 0 
galaxy stellar mass function and cosmic star formation rate 
history, and it has been subsequently shown that it broadly 
reproduces the observed evolving galaxy stellar mass func¬ 
tion out to high redshift (Torrey et al. 2014; Genel et al. 
2014; Sparre et al. 2015). The relatively large volume allows 
for sampling across a range of galaxy environments (Vogels¬ 
berger et al. 2014a) including rare objects (e.g., compact 
massive galaxies, Wellons et al. 2015a) with diverse forma¬ 
tion histories (Sparre et al. 2015), all of which is impor¬ 
tant for the present work. The Illustris simulation contains 
roughly 1820 3 baryon and dark matter particles yielding a 
baryon mass resolution of Mb ar sa 1.3 x 10 6 AFq and a dark 
matter mass resolution of Mdm = 6.3 x 10 6 M© (The number 
of dark matter particles remains exactly fixed at 1820 3 for 
the whole run, but the number of baryonic resolution ele¬ 
ments changes owing to cell (de)refinement). The Plummer- 
equivalent gravitational force softening lengths used in the 
simulation is e = 1.0 /i _1 ckpc for both dark matter and 
baryons until 2 = 1, at which time the baryonic gravita¬ 
tional softening length is capped at a maximum physical 
value of e = 0.5 h~ x pkpc. (The dark matter gravitational 
softening length continues at a fixed comoving size to 2 = 0.) 

Several steps have been taken to post-process the 
Illustris data output to facilitate the present analysis. 
First, the simulation output is run through SUBFIND to 
identify friends-of-friends (FoF) haloes and bound sub¬ 
haloes (Springel et al. 2001; Dolag et al. 2009). Through¬ 
out this paper, we employ the SUBFIND sub-halo catalog 
to identify galaxy populations, including both centrals and 
satellites. Wherever we refer to galaxies or galaxy popula¬ 
tions, we are in detail referring to the self-bound sub-halo 
structures identified by SUBFIND. 

Second, a wide range of physical properties - includ¬ 
ing stellar mass, star formation rate, half-mass radius, etc. 
- of each structure identified with SUBFIND have been tab¬ 
ulated. A catalog of galaxy properties is calculated for each 
galaxy and each redshift independently. Throughout this 
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paper we use stellar masses and dark matter halo masses. 

In both cases, we calculate the stellar (dark matter halo) 
masses as being the total mass of all gravitationally bound 
stellar (dark matter) particles of a given SUBFIND subhalo. 

For this paper, we have additionally calculated the stellar 
velocity dispersion for the galaxy population defined as the 
three dimensional standard deviation of stellar particle ve¬ 
locities calculated within the stellar half-mass radius. « 

The third post-processing step is to link the galaxy cat- ^ 
alogs together in time using merger trees. In this paper we ^ 
adopt the SubLink merger trees as described in Rodriguez- 
Gomez et al. (2015). The merger trees are constructed by 
identifying progenitor/descendant galaxy pairings based on 
overlapping particle compositions identified through particle 
identification numbers. The merger trees facilitate tracking 
of individual galaxies forward and backward in time while 
including in situ growth and contributions from mergers. 

When galaxy mergers occur, the progenitors are segregated 
into a single main branch and secondary progenitors. We de¬ 
fine the main progenitor branch as being the most massive 
branch when summed over the entire formation history un¬ 
til that point (De Lucia & Blaizot 2007; Rodriguez-Gomez 
et al. 2015). Other operational definitions of main progeni¬ 
tor branch are possible (e.g., most massive halo at the previ¬ 
ous snapshot, Springel et al. 2005b) and some of the results 
quoted in this paper depend on this assumption. However, 
we have verified that this choice has a very limited impact 
on our results, with all of our results being qualitatively in¬ 
variant to this choice. 

The full data from the Illustris simulation - including 
all data, post-processing SUBFIND galaxy property catalogs, 
merger trees data, and basic scripts and procedures required 
to reproduce our analysis - have been made publicly avail¬ 
able (Nelson et al. 2015). 1 

3 RESULTS: TRACING GALAXIES VIA 

STELLAR MASS 

3.1 Cumulative Stellar Mass Function 

Perhaps the most relevant aspect of our model for this pa¬ 
per is its ability to reproduce the (cumulative) galaxy stellar 
mass function at many observational epochs. It has also been 
shown that the feedback model employed by Illustris - de¬ 
scribed in detail in Vogelsberger et al. (2013) - is capable of 
producing a galaxy stellar mass function and star formation 
main sequence that broadly matches observations (Torrey 
et al. 2014; Genel et al. 2014; Sparre et al. 2015). This agree¬ 
ment is achieved through a combination of star formation 
driven winds to moderate star formation in low mass galax¬ 
ies, and AGN feedback to regulate the growth of massive 
galaxies. This combination of feedback results in a multi¬ 
epoch galaxy stellar mass function that is similar to mod¬ 
ern semi-analytic models and other hydrodynamical simula¬ 
tions (see Somerville & Dave 2014, for comprehensive review 
plots and discussion). Here, we present fits to the redshift 
evolution of the cumulative galaxy stellar mass function as 
found in our simulations. This fit is important to the analysis 
that we carry out in subsequent sections of the paper. 

1 http://www.illustris-project.org 
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Figure 1. Cumulative stellar mass functions derived from the 
galaxy populations found in Illustris are shown at several red- 
shifts as indicated in the legend. The grey region identifies the 
stellar mass range (vertical strip) and cumulative number density 
range (horizontal strip) that correspond to the Milky Way mass 
objects at redshift z = 0 as defined and discussed in the text. The 
dashed lines shown within indicate the multi-epoch CMF fitting 
functions. The fitting functions nearly overlap with the actual 
CMFs at all redshifts, and so we also show the “error” associated 
with these fits in the panel inset, with the solid blue band indi¬ 
cating 5% errors. The mass evolution of galaxies can be inferred 
from the fitting functions by identifying the mass associated with 
a constant comoving number density at several redshifts (e.g., 
where the grey horizontal band intersects the CMFs). 

Figure 1 shows the cumulative mass function (CMF) at 
several redshifts as realized in the Illustris simulation. We fit 
the simulated cumulative galaxy stellar mass function with 
a power law plus exponential dependence of the form 

N = A exp(—M*) (1) 

where M* = A/*/(10 7 M@). The combined power-law, ex¬ 
ponential form of Equation (1) is adopted to be similar to 
the Schechter (1976) function commonly used to describe 
galaxy stellar mass and luminosity functions. We allow all 
of the fit variables to vary with redshift according to 


A 

= ao + aiz + a, 2 Z 2 

(2) 

a 

— ao + aiz + a2Z 2 

(3) 

0 

= /3 0 + fiiz + /3 2 z 2 

(4) 

7 

= 70 + 712 + 722 2 

(5) 


where 2 is redshift. Adopting this fitting form results in 12 
independent coefficients (all variables on the RHS of equa¬ 
tions 2-5) that are set using an ordinary least squares re¬ 
gression on the CMFs over the redshift range z = 0 to 
z = 6, mass range M* > 10 7 Mq, and number density range 
N > 3 x 10” 5 Mpc _3 . The resulting fits are shown in Fig¬ 
ure 1 as dashed lines, which can be compared against the 
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solid lines that trace the CMFs taken directly from the sim¬ 
ulation. In the inset of Figure 1 we also show Log 10 (AN) = 
Log 10 (iVsi m /iVfit) which gives an impression for the level of 
fit accuracy. Using the power law plus exponential fit de¬ 
scribed above we obtain a fit to the CMF that is valid from 
z — 0 to « = 6 with typical errors of order 1%, which al¬ 
ways remain well below 10%. The best fit coefficients for the 
CMF in Illustris in the redshift range 0 < 2 < 6 are given 
in Table 1. The appropriate limits on this fitting function 
cover the mass range 1O 7 M 0 < M, < 1O 12 M 0 , CMF values 
4> > 3 x 10 -5 Mpc -3 , and redshift range 0 < 2 < 6. How¬ 
ever, one should additionally bear in mind that the baryon 
particle mass in our simulations is ~ 1O 6 AF 0 , and so cau¬ 
tion should be taken when considering the low-mass end of 
the mass function where only 10 — 100 stellar particles are 
included in each galaxy. 

Despite its 12 free terms, Equation (1) is trivial to calcu¬ 
late. Upon evaluation, one can easily identify the cumulative 
number density of galaxies over the full resolved mass range 
in Illustris from 0 < 2 < 6. This is generally useful, includ¬ 
ing for comparisons with observed cumulative stellar mass 
functions - which is outside the scope of this paper. 2 The 
general form allows us to obtain useful fits to a wide variety 
of smoothly varying functions making it possible to use ex¬ 
pressions that take a similar form to Equation (1) at several 
points in this paper including the cumulative velocity disper¬ 
sion function and tracked galaxy number-density evolution. 
We are additionally making available with this paper simple 
python scripts that allow one to evaluate Equation (l). 3 


3.2 Milky Way mass galaxies: constant vs. 
non-constant number density 

We begin our analysis of the evolutionary tracks and evolv¬ 
ing comoving number-densities of Illustris galaxies by con¬ 
sidering the formation history of a population of Milky Way 
mass galaxies. 

We adopt a dehnition for “Milky Way mass galaxies” 
as those galaxies with a redshift 2 = 0 stellar mass in the 
range 4 x 1O 1O M 0 < M* < 5 x 1O 1O M 0 (e.g., McMillan 
2011; Bovy & Rix 2013). This corresponds to ~410 galaxies 
at redshift 2 = 0, including all morphological types, forma¬ 
tion histories, environments, etc. sampled in the simulation 
volume. The vertical grey shaded region in Figure 1 indi¬ 
cates the redshift 2 = 0 mass range adopted for Milky Way 
type galaxies in this section. The corresponding horizontal 
grey shaded region identifies the cumulative number-density 
range that is associated with the redshift 2 = 0 Milky Way 
galaxy mass range. We note that in what follows, our results 
are not very sensitive to this specific choice of initial mass 
or number-density range. 


Table 1. Best-fit parameters to the redshift-dependent CMF 
presented in Equation (1). 



i = 0 

2—1 

i = 2 


-2.893811 

0.082199 

-0.123157 

OLi 

-0.625598 

0.086216 

-0.049033 

Pi 

-0.038895 

0.025419 

-0.007130 

7 i 

11.523852 

-0.187102 

0.021022 


Table 2. Best-fit parameters to the backward-tracked number- 
density evolution. The below parameters along with Equation (1) 
describe the median number density of a redshift 2 = 0 selected 
galaxy population at some higher redshift out to z = 3. 



i = 0 

2—1 

i = 2 


-2.643961 

-0.299579 

-0.037861 

OLi 

-0.526844 

-0.138136 

-0.032246 

Pi 

-0.026482 

-0.016006 

-0.005645 

a 

11.339278 

0.391025 

-0.015732 


Table 3. Best-fit parameters to the forward-tracked number- 
density evolution starting from 2 = 1. The fit is only valid from 
2 = 1 to 2 = 0 . 



i = 0 

2—1 

i = 2 


-3.640099 

1.036010 

-0.311620 

OLi 

-0.797215 

0.292991 

-0.077466 

Pi 

-0.048383 

0.033631 

-0.004868 

a 

11.827591 

-0.665051 

0.184837 


Table 4. Best-fit parameters to the forward-tracked number- 
density evolution starting from 2 = 2. The fit is only valid from 
2 = 2 to 2 = 0. 



i = 0 

2—1 

i = 2 

CLi 

-3.991685 

0.348706 

0.050966 

OLi 

-0.793825 

0.019329 

0.036924 

Pi 

-0.038297 

-0.000232 

0.006948 

a 

11.828429 

-0.280904 

-0.030598 


Table 5. Best-fit parameters to the forward-tracked number- 
density evolution starting from 2 = 3. The fit is only valid from 
2 = 3 to 2 = 0. 



i = 0 

2—1 

i = 2 

Oi 

-4.353012 

-0.438934 

0.229633 

OLi 

-0.792804 

-0.191351 

0.066223 

Pi 

-0.029089 

-0.020006 

0.007395 

7 i 

11.801074 

0.197901 

-0.149748 


2 The CMF is most valuable to the present paper where we focus 
on the evolution of galaxies in number density space. We addition¬ 
ally provide fits of the same form to the differential stellar mass 
function in Appendix A which is the more commonly adopted 
form for examining the galaxy stellar mass function, 
https : //github. com/ptorrey/torrey_cmf 


Table 6. Best-fit parameters to the redshift-dependent cumula¬ 
tive velocity dispersion function presented in Equation (6). 


i = 0 

2 — 1 2 = 2 

CLi 

7.391498 

5.729400 -1.120552 

OLi 

-6.863393 

-5.273271 1.104114 

Pi 

2.852083 

1.255696 -0.286638 

7 i 

0.067032 

-0.048683 0.007648 


© 0000 RAS, MNRAS 000 , 000-000 
































6 P. Torrey et al. 

3.2.1 Examples of Individual Galaxy Evolutionary Tracks 

Figure 2 shows synthetic stellar light images (Torrey et al. 
2015; Snyder et al. 2015) in the SDSS-g, -r, -i bands made 
for a selection of 5 galaxies from the Milky Way mass sample 
at 8 redsliifts during their formation (as labeled within the 
top panel of each column). The progenitors of these 5 sys¬ 
tems have been determined directly from the galaxy merger 
tree which finds the progenitor(s) of any galaxy based on the 
particle ID composition of each galaxy. If multiple progen¬ 
itor galaxies exist while tracing galaxies backward in time 
we always select the “main progenitor”, defined as the pro¬ 
genitor with the most massive history (De Lucia & Blaizot 
2007; Rodriguez-Gomez et al. 2015). We have ordered the 
galaxies by their redshift 2 = 3 progenitor masses, with in¬ 
dividual mass and number density evolution tracks shown 
in the bottom panels. 

We find that three of the systems (the top three) evolve 
without significant influence from mergers, with relatively 
smooth mass growth (and number density tracks). These 
three systems preserve their rank order but diverge in their 
stellar mass from each other with time. By redshift z = 3, 
the top (blue) and middle (green) galaxies have stellar 
masses that are different by an order of magnitude. The 
bottom two examples were selected to highlight systems that 
undergo mergers and change their rank order significantly 
with time. For example, the top (blue) and fourth (yellow) 
galaxies have nearly identical mass evolution tracks from 
redshift z = 0 out to redshift a = 1.5. At that time the 
fourth system can be identified through the postage stamp 
images to be undergoing a significant merger event, after 
which the blue and yellow mass growth tracks diverge. De¬ 
spite their nearly identical mass evolution out to redshift 
1.5, these systems are offset by roughly an order of mag¬ 
nitude in their stellar mass by redshift z = 3. A. similar 
qualitative story holds for the bottom (red) system, which 
follows a relatively median mass growth track out to red¬ 
shift z = 0.7, but quickly becomes the least massive of these 
galaxy progenitors thereafter. 

The mass growth tracks shown in the left bottom panel 
of Figure 2 directly translate via the CMF into the number 
density evolution tracks shown in the bottom right panel. 
We find that the dark and light blue galaxies which were 
already massive systems at redshift 2 = 3 and followed mild 
growth paths thereafter are close to remaining on constant 
comoving number density evolution tracks. In contrast, the 
red system which grew rapidly since redshift 2 = 3 has an 
evolution in number density of nearly an order of magnitude. 
Figure 2 highlights the diversity in individual growth paths 
that occur at a fixed 2 = 0 stellar mass. Using the full sample 
of Milky Way mass galaxies in the simulation, we can further 
consider the median growth tracks and the dispersion about 
those tracks for this galaxy population. 

3.2.2 Population Evolutionary Tracks 

We next consider the mass evolution of the full Milky Way 
mass selected galaxy population using now two complemen¬ 
tary methods. First, armed solely with an evolving set of cu¬ 
mulative stellar mass functions, we can identify the galaxy 
mass associated with a specific number density at any red¬ 
shift. If we assume that progenitor/descendant galaxy pop¬ 


ulations can be matched between different epochs at a con¬ 
stant comoving number density, we can infer the mass of 
Milky Way progenitor galaxies at higher redshifts by con¬ 
sidering where the horizontal grey strip overlaps with the 
CMF at those redshifts as shown in Figure 1. This is the 
method commonly adopted when working with multi-epoch 
extragalactic observational data. Here, this is achieved by in¬ 
verting Equation (1) numerically using a Newton-Raphson 
root finding algorithm to solve for M» = A4,(z,N) given 
some constant choice of the galaxies’ cumulative number 
density, N as a, function of redshift 2 . 

Second, the stellar mass evolution of the simulated 
galaxies can be measured directly from the galaxy merger 
tree, as was demonstrated in the previous subsection. The re¬ 
sults from both mass tracking methods are presented in the 
left panel of Figure 3 with the red line indicating the median 
stellar mass evolution from number-density selection and the 
blue solid line indicating the median stellar mass evolution 
from merger tree analysis. The blue shaded regions show the 
range of stellar masses present from the tracked galaxy pop¬ 
ulation, as noted in legend. We note that, by definition, the 
blue line overlaps identically with the red line at the redshift 
where the galaxy populations is selected. 

Two conclusions can be drawn from the left panel of 
Figure 3: (1) The median mass evolution tracks based on the 
comoving number-density selection and direct merger tree 
tracking are offset from one another and (2) The merger tree 
tracking method possess significant scatter that grows as a 
function of redshift, indicating that there is a large diversity 
in the way galaxies assemble their mass, even when these 
are selected from a relatively narrow mass bin at 2 = 0. 

The stellar mass evolution for the explicitly tracked 
populations is more rapid than the inferred mass evolution 
from number-density selections. By redshift 2 = 2(3) there 
is a factor of ~ 2(4) stellar mass offset between the two 
mass tracks. The small scatter in the number-density se¬ 
lected mass evolution directly results from the assumption 
that galaxies remain rank ordered in their mass at all times. 
In contrast, although all galaxies in the explicitly tracked 
galaxy populations are initially clustered in their stellar 
masses, the range of stellar mass values disperses with time 
as galaxies experience variable growth rates and stochastic 
evolution. 

Also shown in the left panel of Figure 3 is a black 
dashed line corresponding to the Milky Way mass evolu¬ 
tion derived from 3D-HST and CANDELS survey data us¬ 
ing the constant comoving number-density selection in van 
Dokkum et al. (2013, see their equation 1). A similar re¬ 
sult was presented in Patel et al. (2013). We find that - 
though there are some differences in their detailed shape - 
the simulated and observational comoving number-density 
mass evolutions are never separated by more than ~ 20% 
over the redshift range 0 < 2 < 3. This is an indication 
that our simulated cumulative stellar mass function evolves 
quite similarly to the observed cumulative stellar mass func¬ 
tion. Any offset between the van Dokkum et al. (2013) line 
and the simulation red curve are driven by inconsistencies 
between the two sets of mass functions. However, more im¬ 
portantly, both of these mass evolution trajectories are offset 
from the explicitly tracked Milky Way mass evolution - by 
a factor of ~ 2 by redshift 2 = 2. A factor of ~ 2 median 
offset is not very severe (in agreement with the conclusions 
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Figure 2. Mock stellar light images are shown for 5 galaxies taken from the Milky Way mass galaxy sample as described in the main 
text. Each system is shown at 8 redshifts from z = 0 to z = 3 to highlight the variety of formation histories that exist for galaxies 
with similar initial stellar masses. The rows are ordered by increasing 2 = 3 progenitor mass. In some cases (the bottom two rows in 
particular) we find merger events contribute significantly to the growth of these systems. The bottom two panels show the mass (left) 
and cumulative number density (right) evolution for the 5 selected galaxies, with the color of the line corresponding to the border color 
for each image. 


of Leja et al. 2013) and is comparable to other uncertainties 
in stellar mass measurements (e.g., initial mass function un¬ 
certainties, age-dust degeneracies, weakly constrained star 
formation histories, etc.). At the low mass end of the z = 2 
progenitor distribution, however, we find that roughly one 
third of explicitly tracked systems have masses that are off¬ 
set by more than an order of magnitude from the constant 
comoving number-density mass trajectory. 

The origin of the offset between the mass evolution 
tracks is simple: galaxies do not remain exactly rank or¬ 
dered, nor is galaxy number (density) a conserved quantity. 
The right panel of Figure 3 shows the distribution of cumula¬ 
tive number density for the initially selected Milky Way type 


galaxies as they are traced back in redshift. The dispersion 
of the galaxies’ cumulative number densities with redshift 
is indicated through the dark blue bands, as noted in the 
legend. The mass and number-density evolution shown in 
the left and right panels of Figure 3 are exactly interchange¬ 
able as long as we are using the cumulative stellar mass 
function to assign rank order. Thus, we reach the same two 
conclusions from the right panel as we did from the left: 
that the median number density of a tracked galaxy pop¬ 
ulation evolves with time, and that the initially clustered 
population of galaxies disperses with time such that no sin¬ 
gle comoving number-density selection can fully recover the 
initial galaxy sample. Correcting the median offset is fairly 
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Figure 3. (Left) The mass evolution of Illustris Milky Way mass galaxies is shown as a function of redshift for two methods used to 
trace galaxy mass growth with time. The wide blue band indicates the mass distribution for a population of galaxies traced backward 
in time explicitly through the merger tree. The red band indicates the inferred mass evolution by assuming constant comoving number 
density and applying the CMF fitting functions presented in Equation (1) with coefficients provided in Table 1. The black dashed line 
indicates the observational inferred Milky Way mass evolution from van Dokkum et al. (2013) using a constant comoving number-density 
assumption. (Right) The number-density evolution of Milky Way mass galaxies is shown as a function of redshift. The wide blue band 
indicates the distribution of number densities that a tracked galaxy population has when traced backward in time through the merger 
tree. We note that there is both a median evolution with redshift and significant scatter, both of which are discussed in the text. 


straightforward, and we provide a clear procedure for how 
to do so in the next subsection. 


3.3 Fits to the non-constant number density 
evolution across galaxy masses 

3.3.1 Tracing Galaxies Backward in Time 

As described above, Equation (1) provides a redshift- 
dependent fit to the number density as a function of mass 
and redshift. When we performed the regression to deter¬ 
mine the best fit parameters in Section 3.1, we specified a 
set of N, M*, and z points based on data from cumulative 
stellar mass functions (i.e., the N and M* pairings as shown 
in Figure 1) at several redshifts. This is sensible because 
this would be the only information that an observer would 
have on multi-epoch galaxy populations. However, since the 
cumulative stellar mass functions are built independently at 
each redshift, no information regarding the number-density 
evolution of individual galaxies or populations of galaxies is 
retained using this approach. As a result we find that when 
we track a population of galaxies explicitly in time they 
follow a non-constant number-density evolution in time as 
shown in Figure 3. 

Using information on the mass and number-density evo¬ 
lution as a function of redshift we can perform a similar 
regression analysis using Equation (1). However, here we 
want to find the cumulative number-density evolution track 
that best describes the actual tracks taken by galaxies in 


our simulation. To achieve this we take N = N(z) to be 
the number density of each individual galaxy as it evolves in 
time from an initial z = 0 mass of M*,o- The only difference 
between the fit that we perform in this section and what we 
did in the previous subsection is the way that the ( N,M,,z) 
pairings are constructed. Here, rather than using the pair¬ 
ings from the CMF, which are constructed independently at 
each redshift, we use (N pairings constructed from 
the merger tree by tracing each galaxy with redshift z = 0 
mass, backwards in time to obtain the number density 

evolution for every galaxy N = N(z). 

The derived parameters from this fitting procedure are 
given in Table 2. Using this fit we are able to infer the ex¬ 
pected median number density a galaxy population will have 
at some redshift z given its initial z = 0 mass, M* : o- Using 
this fit in conjunction with the tabulated CMF presented in 
Section 3.1 we can then infer the average mass associated 
with this galaxy population at other redshifts. 

We demonstrate the derived mass and number-density 
evolution in the left and right panels of Figure 4, respec¬ 
tively. Figure 4 shows the mass and number-density evolu¬ 
tion of a set of four different galaxy populations where the 
colored bands are constructed from tracking galaxies along 
their main progenitor branch, as described above. All of the 
same behavior that was present for the Milky Way mass 
bin inspected above is also seen for the other galaxy mass 
bins considered here. The black solid lines in each panel 
indicate the mass and number density evolution that is ob¬ 
tained from the constant comoving number density assump- 
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Figure 4. (Left) The mass evolution of z = 0 galaxies tracked backward in time is shown as a function of redshift for three methods used 
to trace the stellar mass growth. As in Figure 3, the wide colored bands indicate the (20, 40, 60, 80)% distribution of galaxies as tracked 
backward in time. (Right) The corresponding number-density evolution is shown for the backward-tracked galaxy populations. There is 
significant scatter in the number-density evolution which becomes worse for the most rare (lowest number-density) bins. In both panels, 
black dashed lines indicate the results of the fit to the median evolution in number density, whose parameters appear in Table 2. In the 
left panel, the inferred mass evolution is obtained by using the CMF to convert the best-fit number density to stellar mass. Similarly, in 
both panels the solid black lines indicate the mass and number density evolution tracks following a constant comoving number density 
assumption. 


tion. The black dashed lines in each panel indicate the mass 
and number-density evolution that is obtained from the non¬ 
constant comoving number density fit. 

By construction, we find that the non-constant number- 
density fit follows the appropriate average trend. We obtain 
the mass evolution shown in the black dashed line left panel 
of Figure 4 by converting the evolving number density into 
a mass based on the tabulated CMF coefficients given in Ta¬ 
ble 1. This procedure can be replicated with observational 
data. The fit parameters from Table 2 are well suited to 
describe the median number-density evolution over the re¬ 
solved mass range M,{z = 0) > 10 9 Mq and redshift range 
0 < z < 3. The fit parameters given in Table 2 can be applied 
to identify the progenitor galaxies that properly follow the 
median mass evolution of an initially selected galaxy pop¬ 
ulation (e.g., Behroozi et al. 2013; Marchesini et al. 2014). 
However, we emphasize that while this fit does describe the 
median number-density evolution, it does not capture the 
scattered growth rates. 

3.3.2 Tracing Galaxies Forward in Time 

An important caveat for the fit parameters presented in Ta¬ 
ble 2 is that they were obtained by identifying galaxies at 
z = 0 and tracing their mass/number density backwards in 
time, and therefore only apply in that direction. This in¬ 
forms us of the mass and number-density evolution of the 
main progenitor galaxies of the galaxy population that is 
present at redshift 2 = 0. However, this analysis is not an in¬ 


herently reversible procedure owning primarily to asymmet¬ 
ric scattered growth rates when examining the forward and 
backward evolutionary paths of galaxy populations. There¬ 
fore, while the results of the previous section can be used 
to identify the past mass or number-density evolution of 
present day selected galaxy populations, we cannot necessar¬ 
ily use the results of the previous subsection to identify the 
present day counterparts to an observed high-redshift galaxy 
population. Instead, to identify the present-day descendants 
of a high-redshift galaxy population, the analogous galaxy 
populations need to be traced forward in time. 

To perform this analysis, we select a galaxy population 
at some non-zero redshift and use the merger trees to follow 
that galaxy population forward in time. Galaxies that merge 
with more massive systems are followed until the merger 
event, after which we assume that the galaxy is no longer 
observable - and so it is not included in the fitting or analysis 
beyond that point. Neglecting these branches entirely does 
not significantly change the results, but the consumption of 
galaxies via mergers can lead to a significant reduction in the 
number of galaxies that are available for forward tracking. 
We find a similar result for the “survival fraction” of galaxies 
to 2 = 0 as has been shown in previous work (e.g., Leja et al. 
2013; Mundy et al. 2015). Specifically, the survival fraction 
is not a very steep function of initial stellar mass, with only 
a weak trend where more massive systems are more likely 
to survive. For a galaxy population selected at redshift 2 = 
2, roughly two-thirds of those galaxies can be expected to 
have redshift 2 = 0 counterparts with the rest having been 
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Figure 5. Same as Figure 4, but with galaxy populations that are tracked forward in time after an initial selection at redshift z = 3. 
(Fit parameters appear in Table 5.) Contrasting this figure with Figure 4 demonstrates the difference between tracking galaxy progenitor 
and galaxy descendant populations. Whereas there is significant median evolution in the number density of backward tracked galaxy 
populations, we find that the overall number density evolution is somewhat weaker for forward tracked galaxy populations. This results 
in the constant comoving number density inferred mass evolution (solid black lines) more closely approximating the median merger tree 
tracked mass evolution (solid colored lines). 


consumed by some larger system. For a z = 3 selected galaxy 
population, the survival fraction drops to roughly half. This 
can be contrasted with the expectation that 100% of redshift 
2 = 0 selected galaxies have meaningful high-redshift main 
progenitors, which only requires assuming that the employed 
simulation has sufficient resolution to continue to track their 
formation backward in time. 

Figure 5 shows the mass and number-density evolution 
of a population of galaxies selected at redshift 2 = 3 and 
tracked forward in time. The mass evolution shown in Fig¬ 
ure 5 is qualitatively consistent with our basic expectations: 
the median galaxy mass increases with time along with an 
increase in the dispersion of the individual galaxy mass dis¬ 
tribution. We perform a regression on the ( N , M* jZ , z) pair¬ 
ings to determine the coefficients of Equation (1), where we 
take M,, z to be the stellar mass of each galaxy at some 
initial redshift 2 (in place of A/*,o from the previous subsec¬ 
tion), and N = N(z) is the mass ranked cumulative number 
density of each galaxy when traced forward in time. We indi¬ 
cate the fit mass and number-density evolution tracks with 
black dashed lines in both panels of Figure 5 and the best fit 
coefficients are given in Tables 3-5. We find that the lower 
three mass bins are tracked very well in time using this fit. 
The highest mass bin shows some significant deviation from 
the fit which is mostly a consequence of the low number of 
galaxies in this bin - which decreases as it moves to lower 
redshift. 

Interestingly, we find that the qualitative behaviors of 
the number-density evolution for the tracked galaxy pop¬ 
ulations are different when tracked forward and backward 
in time. When tracked backward in time (Figure 4), the 


number density of the tracked galaxy population steadily 
increases. However, when tracked forward in time the me¬ 
dian number density of the tracked galaxy population re¬ 
mains much more constant. For comparison, the solid black 
lines in Figure 5 indicate evolutionary tracks of constant co¬ 
moving number density. We find that the median mass and 
number density evolution for one of the bins - the red band, 
which was selected to contain galaxies with stellar masses of 
M* = 1O 1O M0 at redshift 2 = 3- almost identically follows 
the constant comoving number density trajectory. The other 
tracked bins are offset from the constant comoving number 
density track, but remain closer to this constant comoving 
number density track than their backward tracking counter¬ 
parts presented in Figure 4. 

If we compare the resulting mass evolution using the 
backward and forward fits, as shown in Figure 6, we find that 
tracing galaxies forward in time yields a noticeably shallower 
mass evolution. The differences in the for ward/backward 
number-density evolution as well as the offsets in the for¬ 
ward/backward inferred mass growth rates are both primar¬ 
ily driven by the scattered growth rates of galaxies, as we 
discuss in Section 5.3. 


4 RESULTS: TRACING GALAXIES VIA 
STELLAR VELOCITY DISPERSION 

4.1 Cumulative Velocity Dispersion Function 

Velocity dispersion has been advocated as a stable proxy 
for galaxy rank order because of its invariance to growth 


© 0000 RAS, MNRAS 000, 000-000 
























Number Density Evolution 11 



z 


Figure 6. The average mass evolution is shown for several mass 
bins to contrast the results that are obtained from tracking galax¬ 
ies forward and backward in time. Both galaxy merger events and 
asymmetric median galaxy growth rates cause the inferred mass 
evolution to be different at the factor of ~2 level depending on the 
directionality of the galaxy tracking. Note that neither of these 
curves is fundamentally more correct than the other, but rather 
they identify different mass evolution tracks as described in more 
detail in the in Sections 3.3.2 and 5.3. 

via galaxy mergers (e.g., Loeb & Peebles 2003; Bezanson 
et al. 2011). If galaxy growth is driven primarily by mergers 
then the central velocity dispersion will evolve by < 30% 
by redshift z = 3 (Hernquist et al. 1993; Hopkins et al. 
2009). In contrast, if internal changes (e.g., puffing up via 
mass loss from quasars; Fan et al. 2008) dominate over merg¬ 
ers in determining structure evolution of low redshift galaxy 
populations, then the velocity dispersion can increase signif¬ 
icantly. While this point has been examined through semi- 
analytic models in the past (Leja et al. 2013; Mundy et al. 
2015), it has not previously been inspected using numerical 
simulations where the velocity dispersion of galaxies can be 
tracked directly. In parallel with the previous section, here 
we present the multi-epoch cumulative velocity dispersion 
function (CVDF) along with a multi-epoch simple fitting 
function to determine its ability to reliably link galaxy pop¬ 
ulations in time. 

To construct the cumulative velocity dispersion func¬ 
tion, we define the central stellar velocity dispersion, ct*, as 
the three-dimensional standard deviation of the stellar ve¬ 
locities within the stellar half-mass radius. We present the 
CVDF at several redshifts in Figure 7. We employ a fit to 
the CVDF of the form 

N(a ,) = A cr“ + ^ LoBcr * exp(— ct*) ( 6 ) 

where ct* = CT*/10 7 km/sec. We perform a regression anal¬ 
ysis to determine the above coefficients (A, a, /3, and 7 ), 
each of which contains a second order redshift dependence 
as given in Equations 2-5. The derived coefficients are given 



Log((T [km/sec]) 

Figure 7. The cumulative velocity dispersion function (CVDF) 
is shown for several redshifts, as indicated in the legend. In con¬ 
trast to the CMF, the CVDF shows comparatively little evolution 
with redshift after z = 2. Multi-epoch fits given in Equation (6) 
are indicated with dashed lines, with the errors associated with 
these fits indicated in the inset plot. Fit parameters can be found 
in Table 6. 

in Table 6 with the results shown as dashed lines in Fig¬ 
ure 7. The inset axes show the error associated with the 
multi-epoch fit, which is accurate at the < 10 percent level 
between redshifts 0 < a < 6 over the velocity disper¬ 
sion range Log(cr*km/sec) > 1.8 and number density range 
N > 3 x 10 _ 3 Mpc“ 3 . 

We find that there is relatively limited evolution in the 
CVDF from z = 2 to « = 0. This low level of redshift evo¬ 
lution was not present for the CMF for any mass scale. We 
do find that there is evolution in the CVDF beyond redshift 
z = 2 which can be quite significant at all velocity dispersion 
values. 


4.2 Evolutionary Tracks in Velocity Dispersion 

The limited evolution in the CVDF is an intriguing feature 
for comoving number-density analysis. If the CVDF is as¬ 
sembled at early times, then perhaps the central velocity 
dispersion evolution is restricted in time. This would hap¬ 
pen in a scenario where galaxies attain their central velocity 
dispersion at early times without significant evolution there¬ 
after, even in the presence of mass growth. This is expected 
for massive quenched galaxies that assemble at early times 
and retain their internal stellar structure while “puffing up” 
at late times from minor merger events (Naab et al. 2007). 
This scenario is less likely to apply to star forming galaxies. 

We test the velocity dispersion rank ordering by trac¬ 
ing several galaxy populations back in time, similar to what 
was done in the previous section. The results are shown in 
Figure 8 . We could select the exact same systems used in 
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Figure 8. Analogous to Figure 4, but where we are now selecting galaxies in bins of number density according to stellar velocity 
dispersion rather than stellar mass. The initial number density ranges sampled in each bin are chosen to match what was used in 
Figure 4. The legend indicates the stellar mass that corresponds to each number density bin. We find that galaxies do not show very 
significant median evolution in their velocity dispersion when traced out to redshift z = 3. However, we find there is significant spread 
associated both in terms of the velocity dispersion and number-density distribution for this galaxy population when traced backward in 
time. The black dashed lines correspond to the same number-density evolution fit shown in Figure 4 which was constructed using the 
CMF, not a new fit using the CVDF. In the left panel, the CVDF was used to convert number density to velocity dispersion. Despite 
having been constructed using stellar mass, the evolving number density fit also appropriately follows the velocity dispersion evolution. 


the previous section to trace backwards in time. However, 
because there is not a perfect 1:1 correlation between stellar 
mass and velocity dispersion, using the same mass-selected 
galaxy population would introduce a somewhat larger ini¬ 
tial spread in the velocity dispersion-assigned cumulative 
number density (the introduced extra scatter in number 
density is roughly a factor of 3). Given our goal of under¬ 
standing how cumulative number-density selection methods 
are able to trace galaxy populations in time, we instead se¬ 
lect a galaxy population that uses the same initial number- 
density limits as in the previous section, but use the CVDF 
rather than the CMF to assign number density. For the ve¬ 
locity dispersion bin centered around a number density of 
6 x 10 -3 Mpc -3 , this results in the selection of 466 galaxies, 
130 of which were also in the mass-selected galaxy sample 
in the same number density range. Although this is a some¬ 
what different initial galaxy selection, we can consider the 
evolution of this galaxy population in velocity dispersion 
and directly compare the evolution in number-density space 
against what we found in the previous section. 

The velocity dispersion evolution is shown in the left 
panel of Figure 8. The median velocity dispersion evolu¬ 
tion is shown with the solid colored lines, while the shaded 
colored regions identify the spread in the evolving velocity 
dispersion distribution for the initially selected galaxy pop¬ 
ulation as indicated in the legend. We find that the median 
velocity dispersion does not signihcantly change over this pe¬ 
riod of time for any of the bins. A typical change of 0.1-0.2 
dex from redshift z = 0 out to z = 3 is found. For com¬ 


parison, the solid black line indicates the evolution along an 
assumed constant comoving number density trajectory. 

However, we find that the mild evolution of the velocity 
dispersion does not directly translate to the proper recovery 
of the initial galaxy population when selected via their co¬ 
moving number density. The right panel of Figure 8 shows 
the number-density evolution as assigned from the CVDF 
for this tracked galaxy population. Despite the mild evolu¬ 
tion of the velocity dispersion, the divergence of this galaxy 
population in number-density space with time is still sig¬ 
nificant. There is a median offset in the number density of 
this tracked galaxy population that grows with time and 
the scatter of the initial galaxy population reaches order of 
magnitude or larger levels in comoving number density by 
redshift z = 3. The general trend that we find in the evolv¬ 
ing number-density distribution is very similar in terms of 
median offset and scatter growth to that found when we 
used the CMF to trace galaxies in time. To highlight these 
similarities, the black dashed line in Figure 8 is not a new fit 
from this data, but rather the number-density evolution de¬ 
termined in the previous section using the CMF (i.e., the co¬ 
efficients given in Table 2). We find that the median number- 
density evolution that we derived for the CMF applies very 
well to the CVDF number-density evolution. 

Given the applicability of the CMF number-density 
evolution fit, we can consider the inferred velocity disper¬ 
sion evolution. Specifically, we can calculate N = N(z) us¬ 
ing the fits from the previous section and then determine 
cr* = er*(V(z)) via Equation (6) with the coefficients given 
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in Table 6 . The result is shown as the black dashed line 
in the left panel of Figure 8 , which is in good agreement 
with the merger tree tracked velocity dispersion evolution. 
This is a clear indication that there is a mean evolution in 
the number density of galaxies present at nearly the same 
level regardless of whether we employ the CMF or CVDF 
to Milky Way mass systems. 


5 DISCUSSION 

The method of matching galaxies between different epochs 
observationally based on their number density is both widely 
used and reasonably physically justified. As has been shown 
in the past - and as we have confirmed in this paper - the 
errors that are introduced into the inferred stellar mass evo¬ 
lution of galaxies between redshift z = 2 and the present day 
when one uses a constant number-density selection rather 
than the explicit galaxy merger tree mass evolution are not 
catastrophic (i.e. of order ~ 0.3 — 0.5 dex). By neglecting 
the scattered growth histories of galaxies, one can imme¬ 
diately link high and low redshift observed galaxy popula¬ 
tions. This is one of the principal methods employed to in¬ 
fer galaxy mass buildup, size growth, and morphology evo¬ 
lution in past literature. However, such an approach does 
not properly link the vast majority of progenitor and de¬ 
scendant galaxies (Mundy et al. 2015). Depending on the 
elapsed time and mass/number-density bin size, the true 
recovery rate of progenitor/descendant galaxies using a con¬ 
stant comoving number-density selection can easily be of 
order ~ 10 — 30% (Leja et al. 2013; Mundy et al. 2015). 

A crude link does exist between high- and low-redshift 
galaxies in their comoving number density, but this link 
evolves with time and includes significant intrinsic scat¬ 
ter. In this paper we have presented the explicitly tracked 
number-density evolution of galaxies based on a hydrody- 
namical simulation of galaxy formation. We find a median 
offset associated with the growth history of any galaxy 
population when compared against the constant comoving 
number-density selection methods. The magnitude of this 
offset is not the same when tracking galaxies forward and 
backward in time. We find that tracking galaxies forward 
in time yields median mass and number density evolution 
tracks that evolve in better agreement with the constant 
comoving number density than when systems are tracked 
backward in time. We have provided simple fitting functions 
that describe the median number-density evolution - both 
forward and backward in time - that can be applied to obser¬ 
vational studies straightforwardly. Once we adopt a simple 
formulation for the non-constant comoving number-density 
evolution, we can recover the median mass evolution of our 
explicitly traced galaxy population from the CMF alone. We 
encourage this to enter into future observational analysis as 
has been done in Marchesini et al. (2014). 

We have examined the claim that velocity dispersion 
can act as a more robust property for linking galaxies to¬ 
gether in time. By constructing the CVDF, we were able to 
apply an identical analysis to the evolution of the velocity 
dispersion of our tracked galaxy population. Although the 
CVDF itself shows limited evolution from z = 0 to 2 = 2, 
there is still significant evolution in the number density of in¬ 
dividual galaxies as assigned through the CVDF. We found 


that the median evolution in the number density for a pop¬ 
ulation of explicitly-tracked galaxies behaved nearly iden¬ 
tically to what we found when we used the CMF. This il¬ 
lustrates an important point: it indicates that there is an 
underlying driver of galaxy number-density evolution that 
impacts our results regardless of the physical quantity on 
which we perform our galaxy rank ordering. 

5.1 Dependence on Baryon Physics 

How much do the prescriptions derived for the galaxy co¬ 
moving number density evolution in this paper depend on 
the specific physics implementations that we have employed 
in the Illustris simulation? The answer to this question is 
fairly critical, since we have focused only on the number 
density evolution of galaxies as characterized by their bary- 
onic properties, which are subject to influence from poorly 
constrained and crudely modeled sub-grid prescriptions for 
many physical processes. To address this point we consider 
the number density evolution of the dark matter haloes di¬ 
rectly, since they are relatively insensitive to the baryonic 
models. We select four galaxy populations using the same 
number density criteria that were employed in Figures 4 
and 8 . The result of tracking these four galaxy populations 
backward in time is shown in Figure 9. We find that the 
characteristic evolution of dark matter halo masses is fairly 
different from what was found for the stellar mass evolution 
- especially at the high-mass end. At the high-mass end, 
massive galaxies tend to quench owing to the AGN feed¬ 
back prescriptions that have been implemented in our simu¬ 
lation. This leads to relatively flat late-time growth rates for 
the stellar masses as shown in Figure 4. In contrast, no such 
late-time flattening of the halo mass growth rate is present 
in Figure 9. All haloes continue to grow rapidly until the 
present day. 

We next consider what this means for the number den¬ 
sity evolution of this galaxy population as shown in the right 
panel of Figure 9. The black dashed lines indicate the num¬ 
ber density evolution calculated in Section 3 based on galaxy 
stellar mass and are therefore identical to those presented in 
Figures 4 and 8 . We find that - despite the visible differences 
in the stellar mass and halo mass growth trajectories - the 
median number density evolution is nearly identical regard¬ 
less of whether we use stellar mass, stellar velocity disper¬ 
sion, or dark matter halo mass to trace the number density 
of galaxies in time. This allows us to conclude that the im¬ 
plementation of baryonic physics and feedback processes as 
included in our simulations does not dominate the number 
density evolution of galaxies. Rather, the stochastic growth 
rate of the underlying dark matter halo is the primary driver 
of the galaxy number density evolution that we find in our 
simulations. We have repeated this analysis tracking galax¬ 
ies forward in time, and have arrived at the same conclusion. 
For this reason, we consider the derived median forward and 
backward number density evolution trends presented in this 
paper to be robustly tied to the underlying dark matter halo 
growth rates and to be relatively independent of the specific 
implementation of baryon physics/feedback adopted in our 
simulation. 

Nevertheless, we still caution that some of the galaxy 
properties considered in this paper are subject to influence 
from the adopted physics/feedback prescriptions employed 
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in our simulations. For example, the size-mass relation de¬ 
rived for the Illustris galaxy population is shifted to larger 
sizes for low mass galaxies when compared against obser¬ 
vations. This is likely an indication of a shortcoming in ei¬ 
ther our treatment of the ISM equation of state or feed¬ 
back implementation and could impact the derived veloc¬ 
ity dispersion explored in Section 4. Similarly, the galaxy 
stellar mass function obtained within the Illustris simula¬ 
tion broadly agrees with observations across a wide range 
of redshifts (Torrey et al. 2014; Genel et al. 2014; Sparre 
et al. 2015; Somerville & Dave 2014), but differs in detail. It 
is therefore possible that future generations of simulations 
or semi-analytic models that better match the galaxy mass 
distribution consistent with observations could yield some¬ 
what different median or scattered comoving number density 
evolution rates. Although we feel confident that the non¬ 
constant comoving number density fits prescribed in this 
paper are an improvement over the constant comoving num¬ 
ber density assumption commonly applied in the literature, 
the previously mentioned caveats along with those discussed 
in Nelson et al. (2015) should be kept in mind when apply¬ 
ing the evolving cumulative number density fits presented 
in this paper. 

5.2 Additional Parameter Dependencies 

A wide range of galaxy properties beyond the stellar mass, 
velocity dispersion, and dark matter halo mass are tracked 
in our simulations. We can therefore consider the role that 
several other galaxy parameters may play in predicting the 
scatter seen in the galaxy number-density evolution. For ex¬ 
ample, it is reasonable to suspect that the relative late/early 
formation times of galaxies can be distinguished based on 
galaxy color. Given the abundance of basic information we 
have about an observable galaxy population at some red- 
shift (e.g., 2 = 0 galaxy masses, sizes, star formation rates, 
colors, etc.), how deterministically can we predict an indi¬ 
vidual galaxy’s evolutionary history? We have shown in this 
paper that galaxy populations of similar stellar mass will 
have large scatter in their formation histories, and it is not 
immediately clear to what extent we can differentiate be¬ 
tween galaxies that will grow faster or slower compared to 
their peers of similar initial mass by considering additional 
galaxy properties. 

We adopt the most straightforward method to identify 
additional parameter dependencies that follows the same ap¬ 
proach used throughout this paper. Specifically, we perform 
an ordinary linear regression using the redshift 2 = 0 galaxy 
stellar masses, stellar velocity dispersions, sizes, star forma¬ 
tion rates, and g — r galaxy colors. We adopt the stellar 
half-mass radii as a proxy for galaxy size and the g — r color 
is calculated based on the Bruzual & Chariot (2003) stellar 
photometric catalogs as tabulated in Torrey et al. (2015). 
The fit that we apply takes the general form 

4 

Log 10 (JV) = ^^a,^ (7) 

i j =0 

where i is a summation over galaxy properties (i = 0 is 
stellar mass, i = 1 is stellar velocity dispersion, etc.), j is 
a summation over polynomial expansion order, and Cij = 
co,i,j+ci : ijz+C2,ijz 2 are the redshift-dependent coefficients. 


The regression is performed using the tracked number den¬ 
sity of the galaxy population with redshift N = N(z) as well 
as the redshift 2 = 0 galaxy properties. The regression yields 
the best possible fit to the number-density evolution of the 
galaxy population backward in time based on the proper¬ 
ties that are known at redshift 2 = 0. If there are residual 
correlations driving the scatter seen in the number-density 
evolution in Figures 4 and 8 , then they will be captured 
with this fitting procedure. We note that the fourth order 
polynomial in Equation (7) gives fits to the number den¬ 
sity evolution which are equally good (i.e. errors of order a 
few percent) as those given in Section 3 when stellar mass 
is the only parameter considered, which makes this a fair 
comparison. 

Given the number-density evolution fit in Equation (7), 
we can derive the stellar mass evolution using the tabulated 
CMF given in Section 3. We are then able to quantify the 
reduction in the scatter of this fit by considering the error 
in the resulting stellar mass estimates. Figure 10 shows the 
median and standard deviation of the log ratio of the pre¬ 
dicted mass to the actual mass at several redshifts for the 
Milky Way mass selected galaxies using both the mass only 
fits given in Section 3, as well as the multi-parameter fit 
given in Equation (7). We find that the multi-parameter fits 
provide a median error which is similar to the mass-only fit, 
but that the one-sigma standard deviation in the scatter is 
reduced by ~ 0.1 dex. While the multi-parameter fit is an 
improvement over the “mass only” fit, this amounts only to 
a ~ 20% reduction in the scatter. Even with an accounting of 
the galaxy stellar masses, sizes, star formation rates, colors, 
and stellar velocity dispersions entering into our analysis, 
our improved fit still has a 0.3(0.4) dex standard deviation 
by redshift z = 2(3). 

The lack of significantly reduced scatter indicates that a 
direct and unambiguous linking cannot be achieved between 
high and low redshift populations given the simulated galaxy 
stellar masses, sizes, star formation rates, colors, and stellar 
velocity dispersions alone. We do not rule out the possibil¬ 
ity of being able to deterministically connect high redshift 
and low redshift galaxy populations in a direct progenitor- 
descendant link, but our results indicate that this would re¬ 
quire information beyond the quantities explored here. We 
have performed the same exercise tracing galaxies forward 
in time, and found similar results (i.e. a ~ 20 % reduction 
in scatter). While marginally reduced, the scatter is still a 
significant component of the overall mass evolution. 

We caution again that some of the galaxy properties 
considered in this section are subject to influence from 
the adopted physics/feedback prescriptions employed in our 
simulations. We specifically note that although the sim¬ 
ulated galaxy stellar mass function from Illustris broadly 
agrees with observations, the most massive galaxies continue 
to experience intermittent periods of star formation activity 
at late times that can lead to non-zero SFRs and greenish 
galaxy colors. Both of these may adversely impact our ability 
to decompose mass-matched galaxy populations into late- 
and early-forming subsamples. It will therefore be interest¬ 
ing to reconsider this problem using other currently available 
numerical simulations (e.g., Schaye et al. 2015) which em¬ 
ploy different physical/feedback prescriptions (Crain et al. 
2015) or with future generations of large volume galaxy for¬ 
mation simulations or semi-analytic models. 
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Figure 9. Analogous to Figures 4 and 8, but where we are now selecting galaxies in bins of number density according to their dark 
matter subhalo mass rather than stellar mass or velocity dispersion. As in Figure 8, the black dashed lines are number-density evolution 
tracks constructed using the CMF, and the cumulative dark matter mass function was used to convert number density to dark matter 
mass in the left panel. We find that the dark matter mass growth of these systems looks somewhat different from the stellar mass growth 
owing to the lack of quenching. Regardless of this difference, the number density evolution is nearly identical to what we obtained for 
both the central velocity dispersion and stellar mass number density analysis, and the fit correctly tracks the dark matter mass evolution. 


5.3 Progenitor/Descendant Tracking Asymmetry 

We have found in Section 3.3.2 that tracing galaxies forward 
in time yields distinctly shallower inferred mass growth rates 
than tracing galaxies backward in time. A similar manifes¬ 
tation of this effect is the qualitatively different number- 
density evolution for galaxies as they are traced forward 
and backward in time. Physically, tracing galaxies forward 
and backward captures different processes. When tracing 
galaxies forward in time, a significant fraction of the tracked 
galaxy population can be “lost” owing to merger events 
when the galaxy being tracked is swallowed by a more mas¬ 
sive system. The forward tracks therefore roughly capture 
the median mass evolution of the surviving galaxy popula¬ 
tion. Tracing galaxies backward in time contains no analo¬ 
gous loss of systems owing to mergers. By definition, any 
galaxy which exists in the simulation at redshift z = 0 
is a main branch. The backward tracks therefore roughly 
capture the median mass evolution of the main progenitor 
galaxies. Since these two tracking methods capture different 
physical galaxy populations, it should perhaps not surprise 
us that they yield qualitatively similar but quantitatively 
different mass evolution tracks. However, if we select only 
main branches in both the forward and backward tracking 
analysis, we find that a nearly identical bias still persists 
between the inferred mass evolution in each direction. The 
reason for this is that while the forward tracking does in¬ 
deed suffer from a net reduction of tracked galaxies with 
time, the systems which are lost owing to mergers are more- 
or-less randomly sampled from the initial population (there 


are marginal correlations with environment, but these leave 
a non-detectable signal). 

The main effect that drives the difference in forward- 
backward mass tracking is the asymmetric sampling of 
galaxy scattered growth given the initially selected galaxy 
population. A population of galaxies selected at low red- 
shift will naturally contain some subset of galaxies which 
had anomalously fast growth histories (i.e. which originated 
from much lower masses). Although these anomalously fast 
growth histories only apply to a small fraction of the galaxy 
population, the steep nature of the galaxy stellar mass func¬ 
tion implies a much higher abundance of low-mass galaxies 
that are able to follow these tracks. Therefore, there is a 
conditional probability set by the shape of the high-redshift 
galaxy stellar mass function that tends to sample fast growth 
histories when tracing galaxies backward in time. Track¬ 
ing galaxies forward in time yields no similar conditional 
probability. Instead, the primary source of galaxy dispersed 
growth histories is simply the scatter in the galaxy stellar 
mass function. We leave a more formal exploration of the 
various mechanisms that drive galaxy number-density evo¬ 
lution to a future study. 


6 CONCLUSIONS 

In this paper we have studied the stellar mass, central stellar 
velocity dispersion, dark matter halo mass, and correspond¬ 
ing comoving number-density evolution of galaxies using 
the Illustris hydrodynamical galaxy formation simulation. 
We have compared the evolutionary paths of galaxy pop- 
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Figure 10. The median and la standard deviation for the log 
ratio of the predicted mass to the actual tracked stellar mass (from 
the merger tree) as a function of time for z = 0 Milky Way mass 
galaxies. Red lines indicate the predictions when z = 0 stellar 
mass is the only parameter considered, and blue lines indicate the 
predictions when the redshift 2 = 0 stellar velocity dispersions, 
stellar half mass radii, star formation rates, and g — r galaxy 
colors are also included. The scatter is reduced by only ~ 20% 
when these additional parameters are considered. 

ulations obtained by assuming that galaxies preserve their 
number density in time (the so-called constant number den¬ 
sity ansatz) and by directly tracking the simulated galaxies 
backward and forward in time via the available merger trees. 

Our main conclusions are as follows: 

• We provide a simple tabulated function that gives the 
cumulative stellar mass function (CMF) and cumulative 
stellar velocity dispersion function (CVDF) from 2 = 0 to 
2 — 6 in the Illustris Simulation (Equation 1 and Tables 1 
and 6). This simple function can be used - as we do in this 
paper - to infer the stellar mass growth of galaxies at a 
fixed number density. The cumulative stellar mass function 
found in the Illustris Simulation can be compared against 
observations, and we note that previous studies have pre¬ 
sented such a comparison (Torrey et al. 2014; Genel et al. 
2014) with favorable results. The functional fit provided in 
this paper for the differential and cumulative galaxy stellar 
mass function should help facilitate future comparisons with 
simulated data and semi-analytic results. 

• We trace galaxies forward and backward in time us¬ 
ing merger trees from the Illustris simulations and find that 
galaxy populations do not evolve along constant comoving 
number-density tracks. They fail to do so because of the 
combined influence of galaxy mergers and scattered galaxy 
growth rates. We find that galaxies that are initially similar 
in their stellar mass, dark matter mass, or central stellar 
velocity dispersion diverge with time. 

• We find that the central stellar velocity dispersion 


evolves only mildly with redshift owing to the combined 
effects of mass and size growth. Despite the mild velocity 
dispersion evolution we find that velocity dispersion yields 
a number density evolution that is not improved over that 
found for stellar mass or dark matter mass assigned num¬ 
ber density evolution. In fact, we find that the evolution of 
the number-density distribution of galaxies evolves nearly 
identically regardless of whether one uses stellar mass, dark 
matter mass, or central stellar velocity dispersion to assign 
number density. 

• There is a systematic bias between the median mass 
growth rate inferred from constant comoving number- 
density analysis and merger tree analysis that we capture 
in our simulations. This bias is driven by a systematic evo¬ 
lution in the median number density of a galaxy popula¬ 
tion when traced in time. The median offset in stellar mass 
growth histories is only a factor of ~2(4) when tracing Milky 
Way type galaxies out to redshift 2 = 2(3). However, we em¬ 
phasize that this offset is systematic, and can be corrected 
for by accounting for the median number-density evolution 
of galaxies with time. 

• We provide a simple tabulated function that describes 
the number-density evolution of simulated galaxies both for¬ 
ward and backward in time (Equation 1 with Tables 2-5). 
We encourage the use of this simple form in place of the 
widely applied constant comoving number density. While 
the non-constant comoving number-density evolution does 
not capture the scattered growth rates that are present for 
our simulated galaxy population, it does account for the first 
order offset for the median galaxy mass and number density 
evolution. 

• A fundamental asymmetry exists between progenitor 
and descendant tracking. We find that the mass trajecto¬ 
ries identified by following progenitor and descendant galaxy 
populations in time yield an offset of a factor of a few, which 
is systemically biased toward faster growth rates when trac¬ 
ing galaxies backward in time. This implies that the pro¬ 
genitors of Milky Way or other mass galaxies would in fact 
be on average lower in mass than would be implied from a 
constant comoving number-density analysis. This has direct 
implications for quoted, e.g., Milky Way mass progenitor 
mass evolutionary histories in the literature. 

• The scatter in the mass formation histories for any ini¬ 
tially similar galaxy population is large. We show that the 
simulated progenitors of present day Milky Way mass galax¬ 
ies span at z ~ 2 more than one order of magnitude in stel¬ 
lar masses. We apply a regression including several galactic 
properties beyond stellar mass (size, star formation rate, 
galaxy color, and stellar velocity dispersion) and find that 
the error in the mass/number-density evolution can only be 
improved marginally (by ~ 20 %). 

• We argue that the intrinsic scatter in galaxy growth 
rates implies that one cannot unambiguously identify galaxy 
progenitor/descendant populations between different obser¬ 
vational epochs. 


In light of these conclusions, statistical methods for linking 
progenitor and descendant galaxy populations may be better 
suited for observationally deriving galaxy mass, size, and 
morphology evolution. 


© 0000 RAS, MNRAS 000, 000-000 









Number Density Evolution 17 


ACKNOWLEDGEMENTS 

PT acknowledges support from NASA ATP Grant 
NNX14AH35G. SW is supported by the National Science 
Foundation Graduate Research Fellowship under grant num¬ 
ber DGE1144152. FM acknowledges support from the MIT 
UROP program. RM acknowledges support from the DOE 
CSGF under grant number DE-FG02-97ER25308. AP ac¬ 
knowledges support from the HST grant HST-AR-13897. 
Support for C-PM is provided in part by the Miller Insti¬ 
tute for Basic Research in Science, University of California 
Berkeley. VS acknowledges support by the DFG Research 
Centre SFB-881 The Milky Way System through project 
Al, and by the European Research Council under ERC-StG 
EXAGAL-308037. LH acknowledges support from NASA 
grant NNX12AC67G and NSF grant AST-1312095. 

The Illustris-1 simulation was run on the CURIE su¬ 
percomputer at CEA/France as part of PRACE project 
RA0844, and the SuperMUC computer at the Leibniz Com¬ 
puting Centre, Germany, as part of GCS-project pr85je. 
The further simulations were run on the Harvard Odyssey 
and CfA/ITC clusters, the Ranger and Stampede supercom¬ 
puters at the Texas Advanced Computing Center through 
XSEDE, and the Kraken supercomputer at Oak Rridge Na¬ 
tional Laboratory through XSEDE. The analysis presented 
in this paper was conducted on the joint MIT-Harvard com¬ 
puting cluster supported by MKI and FAS. 


REFERENCES 

Aragon-Salamanca, A., Ellis, R. S., Couch, W. J., & Carter, 
D. 1993, MNRAS, 262, 764 

Barro, G., Faber, S. M., Perez-Gonzalez, P. G., et al. 2013, 
ApJ, 765, 104 

Behroozi, P. S., Marchesini, D., Wechsler, R. H., et al. 2013, 
ApJ, 777, L10 

Belli, S., Newman, A. B., & Ellis, R. S. 2014, ApJ, 783, 117 

Bennett, C. L., Halpern, M., Hinshaw, G., et al. 2003, 
ApJS, 148, 1 

Bezanson, R., van Dokkum, P. G., Franx, M., et al. 2011, 
ApJ, 737, L31 

Bovy, .J., & Rix, H.-W. 2013, ApJ, 779, 115 

Brammer, G. B., Whitaker, K. E., van Dokkum, P. G., 
et al. 2011, ApJ, 739, 24 

Bruzual, G., & Chariot, S. 2003, MNRAS, 344, 1000 

Butcher, H., & Oemler, .Jr., A. 1984, ApJ, 285, 426 

Carollo, C. M., Bschorr, T. J., Renzini, A., et al. 2013, ApJ, 
773, 112 

Crain, R. A., Schaye, J., Bower, R. G., et al. 2015, MNRAS, 
450, 1937 

Daddi, E., Renzini, A., Pirzkal, N., et al. 2005, ApJ, 626, 
680 

De Lucia, G., & Blaizot, .J. 2007, MNRAS, 375, 2 

Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 
433, 604 

Dolag, K., Borgani, S., Murante, G., & Springel, V. 2009, 
MNRAS, 399, 497 

Fan, L., Lapi, A., De Zotti, G., & Danese, L. 2008, ApJ, 
689, L101 

Genel, S., Vogelsberger, M., Springel, V., et al. 2014, MN¬ 
RAS, 445, 175 


Guo, Q., White, S., Boylan-Kolchin, M., et al. 2011, MN¬ 
RAS, 413, 101 

Hernquist, L., Spergel, D. N., & Heyl, J. S. 1993, ApJ, 416, 
415 

Hopkins, P. F., Hernquist, L., Cox, T. J., Keres, D., & 
Wuyts, S. 2009, ApJ, 691, 1424 
Jaacks, J., Finkelstein, S. L., & Nagamine, K. 2015, ArXiv 
e-prints 

Katz, N., Weinberg, D. H., & Hernquist, L. 1996, ApJS, 
105, 19 

Keating, S. K., Abraham, R. G., Schiavon, R., et al. 2015, 
ApJ, 798, 26 

Klypin, A. A., Trujillo-Gomez, S., & Primack, J. 2011, ApJ, 
740, 102 

Leja, .J., van Dokkum, P., & Franx, M. 2013, ApJ, 766, 33 
Lidman, C., Suherli, .J., Muzzin, A., et al. 2012, MNRAS, 
427, 550 

Loeb, A., & Peebles, P. J. E. 2003, ApJ, 589, 29 
Marchesini, D., Muzzin, A., Stefanon, M., et al. 2014, ApJ, 
794, 65 

McMillan, P. J. 2011, MNRAS, 414, 2446 
Morishita, T., Ichikawa, T., Noguchi, M., et al. 2015, ApJ, 
805, 34 

Mundy, C. J., Conselice, C. J., & Ownsworth, J. R. 2015, 
ArXiv e-prints: 1504.05583 

Naab, T., Johansson, P. H., Ostriker, J. P., & Efstathiou, 
G. 2007, ApJ, 658, 710 

Nelson, D., Pillepich, A., Genel, S., et al. 2015, ArXiv e- 
prints: 1504.00362 

Papovich, C., Finkelstein, S. L., Ferguson, H. C., Lotz, 
.J. M., & Giavalisco, M. 2011, MNRAS, 412, 1123 
Patel, S. G., van Dokkum, P. G., Franx, M., et al. 2013, 
ApJ, 766, 15 

Rodriguez-Gomez, V., Genel, S., Vogelsberger, M., et al. 
2015, MNRAS, 449, 49 

Saglia, R. P., Sanchez-Blazquez, P., Bonder, R., et al. 2010, 
A&A, 524, A6 

Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 
446, 521 

Schechter, P. 1976, ApJ, 203, 297 

Sijacki, D., Springel, V., Di Matteo, T., & Hernquist, L. 
2007, MNRAS, 380, 877 

Sijacki, D., Vogelsberger, M., Genel, S., et al. 2014, ArXiv 
e-prints: 1408.6842 

Snyder, G. F., Torrey, P., Lotz, J. M., et al. 2015, ArXiv 
e-prints 1502.07747 

Somerville, R. S., & Dave, R. 2014, ArXiv e-prints 
1412.2712 

Sparre, M., Hayward, C. C., Springel, V., et al. 2015, MN¬ 
RAS, 447, 3548 

Springel, V. 2010, MNRAS, 401, 791 
Springel, V., Di Matteo, T., & Hernquist, L. 2005a, MN¬ 
RAS, 361, 776 

Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289 
Springel, V., White, M., & Hernquist, L. 2001, ApJ, 549, 
681 

Springel, V., White, S. D. M., Jenkins, A., et al. 2005b, 
Nature, 435, 629 

Torrey, P., Vogelsberger, M., Genel, S., et al. 2014, MN¬ 
RAS, 438, 1985 

Torrey, P., Snyder, G. F., Vogelsberger, M., et al. 2015, 
MNRAS, 447, 2753 


© 0000 RAS, MNRAS 000, 000-000 



18 P. Torrey et al. 


Trujillo, I., Conselice, C. J., Bundy, K., et al. 2007, MN- 
RAS, 382, 109 

Valentinuzzi, T., Fritz, J., Poggianti, B. M., et al. 2010, 
ApJ, 712, 226 

van Dokkum, P. G., & Franx, M. 1996, MNRAS, 281, 985 

—. 2001, ApJ, 553, 90 

van Dokkum, P. G., Whitaker, K. E., Brammer, G., et al. 
2010, ApJ, 709, 1018 

van Dokkum, P. G., Leja, J., Nelson, E. J., et al. 2013, ApJ, 
771, L35 

Vogelsberger, M., Genel, S., Sijacki, D., et al. 2013, MN¬ 
RAS, 436, 3031 

Vogelsberger, M., Genel, S., Springel, V., et al. 2014a, MN¬ 
RAS, 444, 1518 

—. 2014b, Nature, 509, 177 

Vulcani, B., Bundy, K., Lackner, C., et al. 2014, ApJ, 797, 
62 

Wake, D. A., Nichol, R. C., Eisenstein, D. J., et al. 2006, 
MNRAS, 372, 537 

Wellons, S., Torrey, P., Ma, C.-P., et al. 2015a, MNRAS, 
449, 361 

—. 2015b, in prep 

Wiersma, R. P. C., Schaye, J., Theuns, T., Dalla Vecchia, 
C., & Tornatore, L. 2009, MNRAS, 399, 574 


APPENDIX A: NON-CUMULATIVE GALAXY 
STELLAR MASS FUNCTION 

In Section 3.1 we provided tabulated fits to the cumulative 
galaxy stellar mass function. Although useful for this paper, 
the CMF is less commonly used in the literature compared 
to the (differential) galaxy stellar mass function. Here, we 
provide similar fits to the (differential) galaxy stellar mass 
function from the Illustris simulation that can be used easily 
for comparisons against other simulations or observational 
data sets. We adopt a functional form of 

* - sSs ‘ 4 “ P( - A) <AI) 

where M* = M*/(10 7 Mq) and the fit coefficients are al¬ 
lowed to vary with redshift as described in equations 2-5. 
We identify the best fit coefficients using an ordinary regres¬ 
sion based on the tabulated differential stellar mass func¬ 
tion over the redshift range 0 < 2 < 6 . The galaxy stellar 
mass functions taken directly from the simulations and the 
associated best fits are shown in Figure Al as solid and 
dashed lines respectively. The inset shows the errors asso¬ 
ciated with the fits, which are marginally larger than what 
was found for the CMF. However, the error remains well 
below 10 % for the full resolved redshift, mass, and number 
density. The appropriate limits on this fitting function cover 
the mass range 10 7 AIq < M * < 10 12 Mq, mass function val¬ 
ues 4> > 3 x 10 _ 5 Mpc“‘ ! dex _1 , and redshift range 0 < z < 6 . 
The best fit coefficients can be found in Table Al and a basic 
python script to evaluate the mass functions can be found 
online . 4 



io" nr 
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Figure Al. Galaxy stellar mass functions derived from the 
galaxy populations formed in Illustris are shown at several red- 
shifts as indicated in the legend. The dashed lines shown within 
indicate the galaxy stellar mass function fitting functions. The 
fitting functions approximate the actual galaxy stellar mass func¬ 
tion at all redshifts reasonably well, with the “error” associated 
with these fits in the panel inset. 


Table Al. The best fit parameters to the redshift-dependent 
differential mass function presented in Equation (Al) are given. 



i = 0 

2—1 

i = 2 


-3.082270 

0.091113 

-0.125720 

OLi 

-0.675004 

0.091193 

-0.049466 

ft 

-0.043321 

0.025282 

-0.007046 

7 i 

11.512307 

-0.190260 

0.021313 


4 https://github.com/ptorrey/torrey_cmf 


© 0000 RAS, MNRAS 000, 000-000 



















